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Abstract 

Many recent papers have questioned Irving and Kirkwood's atomistic expression for stress. In 
Irving and Kirkwood's approach both interatomic forces and atomic velocities contribute to stress. 
It is the velocity-dependent part that has been disputed. To help clarify this situation we investigate 
(i) a fluid in a gravitational field and (ii) a steadily rotating solid. For both problems we choose 
conditions where the two stress contributions, potential and kinetic, are significant. The analytic 
force-balance solutions of both these problems agree very well with a smooth-particle interpretation 
of the atomistic Irving-Kirkwood stress tensor. 
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I. INTRODUCTION 



In 2003 Zhou published his lengthy and detailed "New Look at the Atomic Level Virial 
Stress" in the Proceedings of the Royal Society of Londonjl]. He criticized the usual Irving- 
Kirkwood virial expression [2] for the pressure tensor P as a sum of potential and kinetic 
terms. The pressure tensor is the same thing as the comoving corotating momentum flux, 
and is also minus the stress tensor, a = —P. The detailed microscopic Irving-Kirkwood 



approach has been used 



'or more than 50 years in the interpretation of atomistic molecular 



dynamics simulations, j^, ^, 0, Q] Averaged over a homogeneous periodic volume V the Irving- 
Kirkwood expression for the pressure tensor gives: 



-aV = PV = P^V + P^V = + J2^PP^^^ 

i<j 



Here Fij is the force (for simplicity we assume a pairwise-additive potential) exerted on 
Particle i by Particle j, where the vector from j to i is Vij. Particle i, at location with 
mass rrii and momentum pi, obeys Newton's equation of motion. 



Zhou stated that only the tensor force sum, '^{Fr)ij = "^Fijrij, contributes to the stress, 
while the tensor momentum sum, ^{pp/'m)i = YliPiPi/^i)^ does not. 

This idea - including the forces but not the momenta - is not quite so outlandish as it 
seems. In solids, where the longtime average of the particle location is a sensible quantity, 
the virial theorem can be written in a similar tensor form omitting the momenta: 

{PV) = J2{iFR),,) ; R, ^ (r,) . 

i<j 

This form is derived in Section II. C of Reference 4. We use angular brackets here to indicate 
longtime averages. In situations including external forces the tensor force sum must also 
include either {F^^^r)i or {F^^^R)i. 

Subramaniyan and Sun 7|] tested Zhou's ideas with molecular dynamics, heating a model 
atomistic solid subject to a variety of external boundary conditions on the particle coordi- 
nates. Their simulations showed that only the full Irving-Kirkwood pressure tensor, poten- 



tial plus kinetic, was consistent with macroscopic thermodynamics. Liu and Qiu 



recently 
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provided a useful list of references supporting both sides of the question. In addition they 
suggest that the Zhou prescription is correct provided that external fields and rotation are 
not involved. Here we explore those latter two conditions separately and explicitly, showing 
that both [1] an external field (gravity) and [2] condensed-phase rotation can be analyzed 
properly with the Irving-Kirkwood pressure tensor, in a way compatible with macroscopic 
continuum mechanics. This suggests that the original Irving-Kirkwood approach is more 
generally useful than is Zhou's suggested modification of it. 

In order to compute continuous differentiable field variables (density, velocity, energy, 
stress, heat flux, ... ) from atomistic molecular dynamics simulations, for comparison to 
corresponding fields generated by continuum mechanics solutions, we recommend the use 
of "smooth-particle" averages. These correspond to smearing individual particle properties 
over a spatial region of size h, the range of the smooth-particle weighting function, as is 
described in a recent text[9|, summarized in Section II, and applied in Section III. 

Because the derivation of the pressure tensor is familiar, and applies both at and away 
from equilibrium^, 5] we do not repeat that here. Instead, in Sections III and IV, we 
describe and study two specially instructive problems involving gravitational and rotational 
forces. We reserve our conclusions and closing remarks for Section V. 



II. SMOOTH-PARTICLE AVERAGES OF ATOMISTIC PROPERTIES 

Irving and Kirkwood chose to localize particle properties at the particle locations using 
delta functions. Though this is convenient for formal analyses, and even natural for mass 
and momentum, a smoothed or smeared-out particle contribution to potential energy and 
to fluxes often simplifies comparisons with continuum mechanics. The smeared approach 
can provide field variables with two continuous spatial derivatives, as we show below. 

Because "action at a distance" makes the exact location of momentum and energy fluxes 
ambiguous we choose to smear out particle contributions within a spatial region somewhat 
larger in extent than the spacing between particles. We use a local weight function with 
a range h, w{r, h) to convert particle properties to continuum field properties. Consider, 
for example, the density p and the velocity f in a fluid or solid composed of particles 
with individual masses and velocities {mj,fj}. In the smooth-particle approach jol-lioll field 
variables, such as the density and velocity at the point r, are defined as /i-dependent (range- 
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dependent) sums of nearby particle contributions: 

P(^) = - rj\) . 

j 

p[r)v{r) = '^^mjVjw{\r — rj\) . 
j 

The sums include all particles within a distance h of point r. A good feature of this approach 
is that these definitions of density and velocity satisfy the continuity equation, p/p = — V ■ w, 
exactly. Here, as is usual, the dot indicates a comoving time derivative following the motion. 



Lucy was one of the inventors of the smooth-particle approach lOj. For convenience we 



use his form for the weighting function in all of our smooth-particle sums, 
WLucy(kl <h)= Cd{1 - Qx^ + 8x^ - 3x^) ; x = \r\/h. 

This form has two continuous derivatives everywhere. The normalizing prefactor Cd depends 
on the dimensionality D, 

Ci = (5/4/1) ; C2 = (S/vr/i^) ; C3 = {105/16nh^) . 

C is chosen so that the spatial integral of the weight function is unity: 

ph I'h ph 

/ w^{r)2dr = / w'^{r)27rrdr = / w^{r)4:7[r'^dr = 1 . 
Jo Jo Jo 

Lucy's polynomial form is the simplest normalized weight function with a maximum value 
at the origin and two continuous derivatives everywhere. In the following section, where we 
consider the mechanical equilibrium of a two-dimensional fluid in a one- dimensional gravi- 
tational fleld, we compute average values of the pressure tensor using the one-dimensional 
form of Lucy's weight function. 



III. GRAVITATIONAL EQUILIBRATION 

Gravitational equilibration is a problem in which both the potential and kinetic contribu- 
tions to stress can play a role. Where a constant gravitational acceleration acts downward 
in y, the simple force balance equation for mechanical equilibrium is, 

dP/dy = {dP/dp){dp/dy) = —pg . 
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The stationary density profile, p{y), can be found provided that the dependence of pressure 
P on the density p is known. As a simple example problem, chosen to highlight the kinetic 
and potential contributions to the virial, we choose to study the molecular dynamics of an 
atomistic system which closely approximates the isothermal fluid equation of state 

P{p,T) = {py2) + pT ; T^l. 

This equation of state closely corresponds to the virial equation of state for two-dimensional 
particles of unit mass at unit temperature interacting with a "Cusp" potential chosen to 
have a spatial integral of unity: 

</>Cu.p(r <h) = (10/7r/i2)(l -xf ; x= \r\/h ; 



h 

2'Kr(t)cnsv{r)dr = 1 . 



{pl/m) = {pl/m) = kT = 1 . 

We use this cusp interaction for the interparticle forces because the model closely corre- 
sponds to the simple and useful thermodynamic equation of state given above. We choose 
the range of the Cusp pair potential h = 3, so that the deviation of the potential part of the 
pressure tensor from that macroscopic equation of state is of order one percent. 

For periodic two-dimensional systems the virial-theorem expression for the potential part 
of the pressure tensor can be expressed in terms of sums over all A^(A^ — l)/2 pairs of 
interacting particlesj4, 5]. For a hydrostatic fluid, where and are each equal to the 
potential part P* of the hydrostatic pressure P, we have: 

P!.V = 5^(xP,.).<. = PylV = Y.^yFy\<^ = (1/2) ^(P ■ r),<, = P^V . 

For a completely random distribution of particles in the volume V the potential part of the 
pressure is then given by a force integral. The integral can be related to the integral of 
the pair potential using integration by parts. With our particular choice of pair potential 
0, with an integral of unity, and particle mass, unity, the resulting hydrostatic pressure is 
simply half the square of the density: 



P*\/ = (1/2) J2iF ■ 



r)i<j - 
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-[N{N -1)/{AV)] / 27Tr^(j)'dr 



h 

2 a'. 







h 



+ [N{N - 1)/{2V)] / 27Tr(f)dr = N{N - 1)/{2V) 
Jo 

^ Np/2 — ^ ~ (l/2)p^ . 

A snapshot from an isokinetic (constant kinetic temperature) simulation appears in Fig- 
ure 1. 

For convenience we have chosen a situation in which the potential and kinetic parts of 
the pressure are equally important. At unit temperature {kT = 1) and a density of 2 
(p = Nm/V = N/V = 2), we have 

P* ~ {Ny2V) = pV2 = 2 ; = pkT = 2 . 

We choose the gravitational acceleration g so that the "weight" of a column of unit width 
and containing Uy particles is equal to the maximum pressure, 4, at the maximum density, 
p{y = 0) = 2. In this case the mechanical equilibrium force-balance density and pressure 
profiles are: 

(p + l){dp/dy) = -pg P-2 + ln(p/2) = -gy ; 

POO 

P{y) = P%y)+P^{y)=g p(y)dy . 

We test these analytic results against a molecular dynamics simulation carried out 
isothermally 



at a constant temperature of unity. At and below the bottom y = 
of the column we place Gux boundary particles in an area of Sux (corresponding to the 
maximum density, 2). See Figure 1. We also include a short-ranged repulsive force, 

F'^P{y < 0) = -lOOy^ , 

which is applied to those few moving particles which occasionally penetrate the boundary 
at y = 0. 

With periodic boundaries in x and a repulsive boundary at ?/ = a 9216-particle simula- 
tion gives the typical configuration we showed in Figure 1. The corresponding kinetic and 
potential pressure profiles, averaged vertically with Lucy's one-dimensional weight function, 
are compared to the analytic force-balance profile in Figure 2. Evidently the agreement is 
quite good, and would be qualitatively in error were the kinetic contribution to the pressure 
tensor omitted. 
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FIG. 1: Gravitational isothermal equilibrium at unit temperature for UxUy = 96 x 96 = 9216 
moving particles above 6 x 96 = 576 boundary particles fixed at the bottom of the system. The 
width of the system is = 96. The height is unbounded. The field strength g = ^/uy is chosen 
so that the maximum density matches that of the fixed particles at the bottom: p = N/V = 2 at 
y = 0. This snapshot is typical of a long simulation used to calculate the smooth-particle pressure 
profiles shown in Figure 2. In all of the figures dimensionless (or "reduced") units are used. These 
follow from the definitions of unity for the particle mass, Boltzmann's constant, and the length 
and energy scales in the interparticle forces derived from the cusp potential of Sec. II and the 
Hooke's-law potential of Sec. IV. 

IV. ROTATIONAL EQUILIBRATION 

Next we consider the influence of the kinetic pressure on the mechanical equilibrium 
of a rotating solid. We can use molecular dynamics to determine the thermal (velocity- 
dependent) properties of an isolated rotationless crystal. For this study we have chosen a 
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4f 




FIG. 2: Comparison of the observed and analytic pressure profiles for the^gravitational problem 
shown in Figure 1. From top to bottom the three curves are the total (T or Irving-Kirkwood) , 
kinetic (K), and potential or Zhou) contributions to the pressure profile. These observed 
pressure contributions are calculated as smooth-particle averages. The points correspond to the 
analytic expressions from the isothermal equation of state Pt = + Pr = + P- 

nearest-neighbor Hooke's-Law interaction, 

0Hooke = ^(kl - d)^ . 

with the force constant k, characteristic length d, and particle mass m all set equal to 
unity. To make contact with continuum mechanics we write the stress tensor in terms of 
the displacement vector u and elastic constants A and rj: 

XV ■ul + r][{Vu) + {Vuf] , 

where / is the unit tensor, with I^x = lyy = 1 and I^y = lyx = 0. For the nearest-neighbor 
Hooke's-law crystal the Lame constants are known, 

X = rj = a/S/IGk; . 
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Figure 3 



-30 < ( x,y ) < +30 -30 < ( x,y ) < +30 

FIG. 3: Stationary rotafion snapshots of two 2335-particle Hboke ^hsw crystals. In the rotation- 
less stress-free case all 6828 nearest-neighbor distances are unity. In the steady rotational situations 
shown here, both with an angular frequency uj = 0.01, the tensile strain offsetting the centrifugal 
forces is maximized at the center of the rotating solid. The left view is a cold solid. The right view 
has a temperature kT = 0.01. 

as is also the complete vibrational frequency distribution along with the bulk and surface 
entropies. See Chapter 4 of Reference 5 for details. 

The radial displacement in a rotating disk of radius R, u{r), as well as the corresponding 



stress tensor a are well-known results of linear elastic theory 



111 ]. A derivation for our 



two-dimensional situation is sketched in the Appendix. The results are: 

u{r) = (cjV/18)[5/?2-2r2] ; 

arr = (pcuVl2)[5i?2 - Sr^] ; age = (pcuVl2)[5i?2 - Sr^] . 

The stress components satisfy the radial force-balance equation for a plane-polar-coordinate 
volume element rdrdO rotating at the angular frequency uj: 

+pr = -pro/ = (darr/dr) + {arr - cree)/r . 



9 (March 3, 2009) 



In the comoving and corotating frame, where stress is the negative of the momentum 
flux, rotation provides a centrifugal force per unit mass varying as lJ^ . 

To compare these results from linear elasticity to molecular dynamics simulations, con- 
sider the stationary rotation of a Hooke's-Law lattice. Figure 3 shows two nominally 
stationary states of a 2335-particle solid with an angular velocity oi uj = 0.01. The cold 
crystal is shown at the left. The kinetic temperature of the warm crystal shown on the right 
is kT = 0.01. The 2335-particle crystal is nearly circular. It is the smallest with 36 parti- 



cles equidistant from the origin (at v637 ~ 25.239). Both these rotational problems were 
initialized by thermostating the radial momenta % |5[ while rescaling the angular momenta 
to generate thermally-equilibrated steadily-rotating solid disks. During the first half of each 
run two separate rescaling, or "Gaussian" , thermostats were applied, so as to keep the radial 
temperature and the angular velocity constant. 

Figure 4 illustrates the approximately-quadratic dependence of the maximum tensile 
stress on the rate of rotation for small angular velocities. For comparison with the simulation 
results the linear-elastic stress at the center of a disk with the same mass, Nm = 2335, and 
a series of rotation rates u is also shown. The agreement is correct to four figures as u 0. 

Let us next consider the stresses in a thermally excited rotating crystal, computed accord- 
ing to the virial theorem using Irving and Kirkwood's formulation of the atomistic stresses. 
The Hooke's-law nature of the particle interactions guarantees that our model crystals will 
not melt. But as temperature rises the deformation can become quite large, so that linear 
elastic theory no longer applies. Figure 5 is a typical view of a rotating specimen at a rotation 
rate of u = 0.01 and a kinetic temperature (relative to rigid-body rotation) kT = 0.02. 

The simplest route to the polar-coordinate stress tensor is, first, to calculate the kinetic 
and potential parts of each particle's pressure tensor in Cartesian coordinates: 

iP!.V), = \ Y.(^xF/r),, ; (P* r), = i 5^(xyP/r),, ; (P^^^V), = ^ ^^(yyP/r),, . 
j j j 

In keeping with the Irving-Kirkwood picture, the potential contributions to the pressure 

tensor are divided evenly between pairs {i,j} of interacting particles. The polar-coordinate 

representation for each particle's pressure tensor follows from the Cartesian representation 

by a simple rotation, which can be written as a pair of matrix multiplications: 

(-P^)polar = R ■ (-P^) Cartesian ' -R* ; 
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0.03 

a(cjo) 

0.02 



0.01 - 



O.OO 



0.0000 < co2< 0.0001 



FIG. 4: Angular velocity dependence of the cold-crystal maximum tensile stress on rotation rate. 
The molecular dynamics data, shown here as points, for nearly circular solids of the type shown in 
Figure 3, agree with the linear elastic rcsTilt (shown as a straight line in the figure) for disks to four 
figures as the rotation rate goes to zero. The linear-elastic result is (JraaxV/N = {bN y/^Ji/12'K)u^ . 




Figure 5 illustrates a thermally- excited rotating Hooke's-law crystal. For the figure we 
have chosen the temperature so that the thermal stresses make a significant contribution to 
the pressure tensor. The radial stress vanishes at the disk boundary, while the circumferential 
"hoop" stress remains tensile there in conformity to the predictions of linear elasticity. 

The stresses in two rotating crystals, one cold and one hot, are compared with the the- 
oretical results from elastic theory in Figures 6 and 7. The agreeement is nearly perfect, 
and would be spoiled if the kinetic contributions were not included. In particular, omitting 
the kinetic contribution to the radial stress would be quite inconsistent with the vanishing 
of that stress component at the boundary of the disk. 




Ri 



-\-cos{6i) -l-sin(^j) 
— sin(^j) -I- cos(^j) 



; 6i = arctan(7//a;)j . 
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Figure 5 

30: 




FIG. 5: View of a rotating 2335-particle Hooke's-law crystal at an angular velocity of 0.01 and a 
temperature kT = 0.02. 



0.01 

0.00 

-0.01 - 

-0.02 - 

-0.03 

< r < 30 

FIG. 6: PV in the rotating cold crystal of Figure 3 with oj = 0.01. The theoretical radial and 
circumferential components are shown as lines based on the expressions derived in the Appendix. 
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0.02 



0.01 



CO= 0.01, kT = 0.01, N = 2335 



0.00 - 



-0.03 



-0.01 - 



-0.02 - 




O < r < 30 



FIG. 7: Time-averaged stresses in the warm rotating thermally-excited crystal of Figure 3 with 
u = 0.01 and kT = 0.01. The thermal contributions to (PF)„- and {PV)gg are the points at the 
top. The theoretical expressions for the stress (based on the cold-crystal elastic constant) shown as 
lines in the Figure agree well with the points representing results from molecular dynamics. The 
molecular dynamics results include both the potential and kinetic contributions to the comoving 
corotating stresses. 

V. CONCLUSION 

Both the gravitational and the rotational problems show excellent correspondence be- 
tween conventional continuum mechanics and atomistic mechanics provided that both the 
kinetic and potential parts of the pressure tensor are included in the analysis. Although for 
stationary solids the solely potential form for the virial theorem is correct, the number and 
type of problems which can be studied numerically is greatly enhanced by including Irving 
and Kirkwood's ideas coupled with the smooth-particle averaging introduced by Lucy and 
Monaghan in 1977. For well-defined local properties, both at, and especially away from 
equilibrium, it is essential that these properties be measured in a coordinate frame that 
moves with the material. It is no accident that the fundamental equations of continuum 
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mechanics take their simplest form in the comoving frame. In particular, the pressure (or 
stress) and temperature tensors, as well as the heat flux, only make sense in this frame. 
Stress and pressure cannot depend upon the chosen coordinate system. Hence we must 

choose the "comoving" "corotating" "Lagrangian" frame. In that frame the pressure tensor 
is simply the momentum flux, and has both potential and kinetic contributions, as shown 
clearly in the two problems solved here. 

VI. APPENDIX 

The stationary rotation, at angular velocity a;, of an elastic disk of radius R with equal 
Lame constants X — rj = y^3/16 obeys the force-balance equation in the comoving frame, 

= ^-fxJ^r + {durr/dr) + {Urr — cFee)/^ ■ 

This macroscopic problem corresponds to a microscopic model composed of unit-mass par- 
ticles linked by nearest-neighbor Hooke's-law springs. Both the spring constant and the 
rest length of the springs are taken equal to unity. The displacement responsible for the 
radial strain — (durr/dr) causes a corresponding strain in the circumferential direction, 
^ee — {u/r). The stresses, 

arr = r][3{du/dr) + (u/r)] ; aee = r][{du/dr) + S{u/r)] , 

convert the force-balance to an ordinary differential equation: 

r^{d'^u/dr'^) + r{du/dr) -u^ -a;V/3 , 

with a unique solution such that the radial stress vanishes at R: 

u{r) = {puj'^r / ASri)[bR^ - 2r^] = (cuV/18)[5i?=^ - 2r^] . 

This solution can be used to generate the maximum tensile stress in the disk as well as the 
stress and strain profiles. 

arr = {puj'^/12)[5R^ - Sr^] ; aee = {puj'^/12)[5R^ - 3r^] . 
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